Midterm Exam
Kelly Luis
1) Sampling your system
Each of you has a study system your work in and a question of interest. Give an example of one variable that you would sample in order to get a sense of its variation in nature. Describe, in detail, how you would sample for the population of that variable in order to understand its distribution. Questions to consider include, but are not limited to: Just what is your sample versus your population? What would your sampling design be? Why would you design it that particular way? What are potential confounding influences of both sampling technique and sample design that you need to be careful to avoid? What statistical distribution might the variable take, and why?
Variable: Remote Sensing Reflectance Remote sensing reflectance (1/m) is the spectral distribution of visible light from below the ocean surface.
Just what is your sample versus your population? A “sample” of of remote sensing reflectance would be a measurement from a given time point and a “population” would be measurements from all time points and from every location possible.
What would your sampling design be? Our group primarily measures remote sensing reflectance with a hyperspectral floating spectrometer from the side of a boat. The floating spectrometer takes measurements of the downwelling and upwelling light fleld in the visible range and the ratio between the two is the remote sensing reflectance. Measurements every few nanometers every few seconds. We deploy our instrument over diverse shallow benthic habitats to determine the remote sensing reflectances characteristics of varying benthic habitat.
Why would you design it that particular way? Our approach is designed this way because very few light field measurements are taken over diverse shallow water environments. Creating a library of varying benthic types, depths, water properties, and sky conditions informs ocean color algorithm models used fror calibrating and validating airborne and shipborne ocean color sensors. Thus, our sampling design involves looking at benthic habitat maps, talking with local scientists about where we can find varying types, and then working with boat captains to see if we are able to take measurements in those areas.
What are potential confounding influences of both sampling technique and sample design that you need to be careful to avoid? Taking measurements during foul weather or non-clear sky days will introduce variability into our measurements of the in-water and above water light field. In addition, in averaging over nanometer and seconds intervals will also introduce uncertainty.
What statistical distribution might the variable take, and why? If I’m looking at measurements of remote sensing reflectance over a single bottom type (let’s use white sand as an example), my variable might have a normal distribution. If I was to take measurements over white sand at different locations and at different time points, the reflectance would center around the mean and I think the bulk of the observations would fall within a standard deviation of the mean. I think if it varied significantly folks in my field would classify the bottom as a different type of benthic habitat (i.e. black sand, mixture, etc.). Thinking about this problem makes me want to look at the distribution of sand remote sensing reflectances.
2) Let’s get philosophical
Are you a frequentist, likelihoodist, or Bayesian? Why? Include in your answer why you prefer the inferential tools (e.g. confidence intervals, test statistics, posterior probabilities, etc.) of your chosen worldview and why you do not like the ones of the other one. This includes defining just what those different tools mean! Kelly’s Worldview: Frequentist…for now
My scientific worldview is tainted by the question: “How well does an airborne or satellite radiometric measurement match a shipborne radiometric measurement?” Thus, my field operates under a frequentist mindset where a true value, specifically the shipborne measurement, exists and inferential tools like confidence intervals are widely used to determine how “well” an observation (i.e. airborne or spaceborne measurement) captured the true value. I can definitley see Bayesian and likelihood methods applied to ocean color inversion models. For example, Frouin & Pelletier (2013) adopt a Bayesian framework to evaluate the uncertainties from marine reflectance retrievals from satellite measurements. Priors (i.e. weather, cloud cover, land cover) are especially useful in constraining ocean color inversion models and accounting for the probability of observing marine reflectance for each model parameter in an inversion would help with interpreting uncertainties. In addition, I could see myself being a likelihoodist because maximum likelihood estimates attempt to best estimate the true value of a parameter. Ultimately, I use frequentist methods because operating my lab and field primairly operate under frequentist methods, but depending on where my path takes me I could see myself dipping into likelihood and bayesian methods.
3) Power
We have a lot of aspects of the sample of data that we collect which can alter the power of our linear regressions. 1. Slope 2. Intercept 3. Residual variance 4. Sample Size 5. Range of X values
Choose three of the above properties and demonstrate how they alter power of an F-test from a linear regression using at least three different alpha levels (more if you want!) As a baseline of the parameters, let’s use the information from the seal data:
slope = 0.00237, intercept=115.767, sigma = 5.6805, range of seal ages = 958 to 8353, or, if you prefer, seal ages ∼ N(3730.246, 1293.485). Your call what distribution to use for seal age simulation.
# Set Parameters
slopeSeal <- 0.00237
interceptSeal <- 115.767
sigmaSeal <- 5.6805
minSeal <- 958
maxSeal <- 8353
sizeSeal <- 10 # Arbitrarypval_ftest <- function(samp_size = sizeSeal, slope = slopeSeal, intercept = interceptSeal, sigma = sigmaSeal, min = minSeal, max = maxSeal){
# Generate Normal Distribution for Seal Age (Days)
age.days <- runif(samp_size, min, max)
# Generate Data Frame for Seal Age and Seal Length (cm)
seal_df <- data.frame(intercept = intercept, sigma = sigma, age.days = age.days,
length.cm = rnorm(samp_size, intercept + slope * age.days, sigma))
# Compute Linear Regression
seal_reg <- lm(length.cm~age.days, seal_df)
# Compute ANOVA
seal_anova <- anova(seal_reg)
# Extract F-Test P-Value
seal_pval <- seal_anova$`Pr(>F)`[1]
return(seal_pval)
}pow_ftest <- function(nsims = 100, alpha =0.05, samp_size = sizeSeal, slope = slopeSeal, intercept = interceptSeal, sigma = sigmaSeal, min = minSeal, max = maxSeal){
# Apply seal_pval nsims times
p <- replicate(nsims, pval_ftest(samp_size, slope, intercept, sigma, min, max))
# Calculate the number of p values that are incorrect given that we should be rejecting the null
num_wrong <- sum(p > alpha)
# Return power
power <- 1 - num_wrong/nsims
return(power)
}pow_df <- crossing(alpha_diff = seq(from = 0.01, to = 0.1, by = 0.01),
size_diff = 5:25) %>%
rowwise() %>%
mutate(power = pow_ftest(alpha = alpha_diff,
samp_size = size_diff)) %>%
ungroup()
size_plot <- ggplot(pow_df, aes(x=size_diff, y = power, color = factor(alpha_diff))) +
geom_point() +
geom_line() +
scale_color_manual(values = beyonce_palette(79))
size_plot
pow_df <- crossing(alpha_diff = seq(from = 0.01, to = 0.1, by = 0.01),
intercept_diff = 110:130) %>%
rowwise() %>%
mutate(power = pow_ftest(alpha = alpha_diff,
intercept = intercept_diff)) %>%
ungroup()
intercept_plot <- ggplot(pow_df, aes(x=intercept_diff, y = power, color = factor(alpha_diff))) +
geom_point() +
geom_line() +
scale_color_manual(values = beyonce_palette(79))
intercept_plot
pow_df <- crossing(alpha_diff = seq(from = 0.01, to = 0.1, by = 0.01),
slope_diff = seq(0.001,0.01,by=0.001)) %>%
rowwise() %>%
mutate(power = pow_ftest(alpha = alpha_diff,
slope = slope_diff)) %>%
ungroup
slope_plot <- ggplot(pow_df, aes(x=slope_diff, y = power, color = factor(alpha_diff))) +
geom_point() +
geom_line() +
scale_color_manual(values = beyonce_palette(79))
slope_plot
General Observations for Exercise 3 1) Increasing the sample size increases the power 2) Increasing the intercept has minimal impact on the power 3) Increasing the slope increases the slope 4) Increasing the alpha increases power 5) As you’ve noticed, there’s some noise along these plot lines. I’ve thought about smoothing these lines, but what is the value of smoothing power lines besides looking pretty?
4) Bayes Theorem
I’ve referenced the following figure a few times. I’d like you to demonstrate your understanding of Bayes Theorem by hand showing what the probability of the sun exploding is given the data. Assume that your prior probability that the sun explodes is p(Sun Explodes) = 0.0001. The rest of the information you need is in the cartoon!
sun_explodes <- 0.0001
detector_right <- 35/36
sun_stable <- 1-sun_explodes
detector_wrong <- 1/36
yes <- (sun_explodes*detector_right)+(sun_stable*detector_wrong)
detector_sun_explodes <- (sun_explodes*detector_right)/(sun_explodes*detector_right)+(sun_stable*detector_wrong)
detector_sun_explodes <- (sun_explodes*detector_right)/yes5) Quailing at the Propect of Linear Models
I’d like us to walk through the three different ‘engines’ that we have learned about to fit linear models. To motivate this, we’ll look at Burness et al.’s 2012 study “Post-hatch heat warms adult beaks: irreversible physiological plasticity in Japanese quail http://rspb.royalsocietypublishing.org/content/280/1767/20131436.short the data for which they have made available at Data Dryad at http://datadryad.org/resource/doi:10.5061/dryad.gs661. We’ll be looking at the morphology data.
5.1 Three fits
To begin with, I’d like you to fit the relationship that describes how Tarsus (leg) length predicts upper beak (Culmen) length. Fit this relationship using least squares, likelihood, and Bayesian techniques. For each fit, demonstrate that the necessary assumptions have been met. Note, functions used to fit with likelihood and Bayes may or may not behave well when fed NAs. So look out for those errors.
# Read Data
quail_raw <- read_csv("Morphology data.csv")## Parsed with column specification:
## cols(
## `Bird #` = col_integer(),
## Sex = col_character(),
## `Age (days)` = col_integer(),
## `Exp. Temp. (degree C)` = col_integer(),
## `Mass (g)` = col_double(),
## `Tarsus (mm)` = col_double(),
## `Culmen (mm)` = col_double(),
## `Depth (mm)` = col_double(),
## `Width (mm)` = col_double(),
## NOTES = col_character()
## )# Clean Column Names
colnames(quail_raw) <- gsub("[()]", "", gsub(" ", "_", colnames(quail_raw)))
# Delete Row with NAs
quail <- quail_raw[!is.na(quail_raw$Tarsus_mm * quail_raw$Culmen_mm),]
# Plot Data
quail_plot <- ggplot(data = quail, aes(x = Tarsus_mm, y = Culmen_mm)) +
geom_point() +
stat_smooth(method = lm) +
xlab("Tarsus (mm)") +
ylab("Culmen (mm)") +
theme_bw()
quail_plot
# Linear Regression Fit
ls_fit <- lm(Culmen_mm ~ Tarsus_mm, data = quail)
# Prediction
predFrame <- data.frame(Tarsus_mm = min(quail$Tarsus_mm):max(quail$Tarsus_mm))
predVals <- predict(ls_fit, newdata=predFrame, interval="prediction")
predFrame <- cbind(predFrame, predVals) %>%
rename(Culmen_mm = fit)
# Plot Predicted Values from Least Squares Fit
ls_plot <- quail_plot +
stat_smooth(method="lm", color="blue") +
geom_ribbon(data=predFrame, mapping=aes(x=Tarsus_mm, ymin=lwr, ymax=upr),
fill="grey", alpha=0.6) +
theme_bw()
ls_plot
# Evaluate Assumptions
# Check Cook's Distance for Values Greater than 1
plot(ls_fit, which=c(4,5))
# Test for Normality
shapiro.test(residuals(ls_fit))##
## Shapiro-Wilk normality test
##
## data: residuals(ls_fit)
## W = 0.99384, p-value = 0.003158lh_fit <- glm(Culmen_mm~Tarsus_mm,
family = gaussian(link = "identity"),
data = quail)
lh_plot <- ggplot(data = quail, aes(x = Tarsus_mm, y = Culmen_mm)) +
geom_point() +
stat_smooth(method = "glm", method.args = list(family = gaussian(link="identity")))+
xlab("Tarsus (mm)") +
ylab("Culmen (mm)") +
theme_bw()
lh_plot
# Evaluate Assumptions
# Fitted vs. Observed
lh_fitobs <- ggplot() +
aes(x=predict(lh_fit), y=residuals(lh_fit)) +
geom_point() +
xlab("Fitted") + ylab("Observed") +
theme_bw(base_size=17)
lh_fitobs
# QQ Plot
qqnorm(residuals(lh_fit), cex.lab=1.5)
qqline(residuals(lh_fit))
# Profile Plot
plot(profile(lh_fit))
# Model Comparison
lh_fit_null <- glm(Culmen_mm ~ 1,
family = gaussian(link = "identity"),
data=quail)
anova(lh_fit_null, lh_fit, test = "LRT")# Check Cook's Distance
plot(lh_fit, which= c(4,5))
# Check for Normality
shapiro.test(residuals(lh_fit))##
## Shapiro-Wilk normality test
##
## data: residuals(lh_fit)
## W = 0.99384, p-value = 0.003158# Bayes Setup
options(mc.cores = parallel::detectCores())
set.seed(607)
# Bayes Fit with brm
bayes_fit <- brm(Culmen_mm~Tarsus_mm,
data = quail,
family=gaussian())## Compiling the C++ model## In file included from file6122b2018c1.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Core:531:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122b2018c1.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:2:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/LU:47:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122b2018c1.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:3:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Cholesky:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Jacobi:29:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122b2018c1.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:3:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Cholesky:43:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122b2018c1.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/QR:17:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Householder:27:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122b2018c1.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:5:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/SVD:48:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122b2018c1.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Geometry:58:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122b2018c1.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:7:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Eigenvalues:58:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122b2018c1.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:96:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Sparse:26:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/SparseCore:66:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122b2018c1.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:96:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Sparse:27:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/OrderingMethods:71:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122b2018c1.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:96:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Sparse:29:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/SparseCholesky:43:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122b2018c1.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:96:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Sparse:32:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/SparseQR:35:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122b2018c1.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:96:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Sparse:33:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/IterativeLinearSolvers:46:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122b2018c1.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:44:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp:14:13: warning: unused function 'set_zero_all_adjoints' [-Wunused-function]
## static void set_zero_all_adjoints() {
## ^
## In file included from file6122b2018c1.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:70:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/autocorrelation.hpp:18:8: warning: function 'fft_next_good_size' is not needed and will not be emitted [-Wunneeded-internal-declaration]
## size_t fft_next_good_size(size_t N) {
## ^
## 15 warnings generated.## Start sampling# Inspect chains and posteriors
plot(bayes_fit)
#Inspect rhat
mcmc_rhat(rhat(bayes_fit))
#Inspect Autocorrelation
mcmc_acf(as.data.frame(bayes_fit))
#model assumptions
bayes_fitted <- predict(bayes_fit) %>% as_tibble
bayes_res <- residuals(bayes_fit)%>% as_tibble
qplot(bayes_res$Estimate, bayes_fitted$Estimate)
#fit
pp_check(bayes_fit, type="scatter")## Using 10 posterior samples for ppc type 'scatter' by default.
#normality - Check for outliers
qqnorm(bayes_res$Estimate)
qqline(bayes_res$Estimate)
pp_check(bayes_fit, type="error_hist", binwidth = 2)## Using 10 posterior samples for ppc type 'error_hist' by default.
##match to posterior - we want a nice cluster here
pp_check(bayes_fit, type="stat_2d")## Using all posterior samples for ppc type 'stat_2d' by default.
pp_check(bayes_fit)## Using 10 posterior samples for ppc type 'dens_overlay' by default.
#coefficients
summary(bayes_fit, digits=5)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Culmen_mm ~ Tarsus_mm
## Data: quail (Number of observations: 766)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept -0.10 0.21 -0.51 0.30 3726 1.00
## Tarsus_mm 0.37 0.01 0.36 0.39 3658 1.00
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 1.24 0.03 1.18 1.30 3493 1.00
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).#confidence intervals
posterior_interval(bayes_fit)## 2.5% 97.5%
## b_Intercept -0.5050770 0.3020018
## b_Tarsus_mm 0.3602278 0.3855762
## sigma 1.1794704 1.3022714
## lp__ -1260.0431710 -1255.4884869# Plot Bayes Fit
bayes_chains <- as.data.frame(bayes_fit)
bayes_plot <- quail_plot +
geom_abline(intercept=bayes_chains[,1], slope = bayes_chains[,2], alpha=0.1, color="lightgrey") +
geom_abline(intercept=fixef(bayes_fit)[1], slope = fixef(bayes_fit)[2], color="red") +
geom_point()
bayes_plot
# Visualize Coefficients
bayes_fit %>%
gather_draws(b_Intercept, b_Tarsus_mm, sigma) %>%
ggplot(aes(x = .value, y = .variable)) +
geom_halfeyeh( .width = c(0.8, 0.95)) +
ylab("") +
ggtitle("Posterior Medians with 80% and 95% Credible Intervals")
# Differences by Chain
bayes_fit %>%
gather_draws(b_Tarsus_mm) %>%
ggplot(aes(x = .value, y = .chain, fill = factor(.chain))) +
geom_density_ridges(alpha = 0.5)## Picking joint bandwidth of 0.0015
General Observations for 5.1 1) ls_fit = residuals, shapiro test (p<0.05), and cook’s distance (<1) ensured normality and homeostacity 2) lh_fit = residuals, shapiro test (p<0.05), and cook’s distance (<1) ensured normality and homeostacity 3) bayes_fit = rhat (1), qqplot (no large or systematic deviations), pp_check (nice cluster), adn ensured normality and homeostacity
5.2 Three interpretations
OK, now that we have fits, take a look! Do the coefficients and their associated measures of error in their estimation match? How would we interpret the results from these different analyses differently? Or would we? Note, confint works on lm objects as well.
# Least Squares Fit Coefficients
ls_metrics <- broom::tidy(ls_fit)
ls_metrics# Anova
ls_anova <- broom::tidy(anova(ls_fit))
ls_anova# Confidence Intervals
ls_conf <- confint(ls_fit)
ls_conf## 2.5 % 97.5 %
## (Intercept) -0.5216505 0.3242363
## Tarsus_mm 0.3598809 0.3859727# Likelihood Fit Coefficients for Intercept Estimate
lh_metrics <- broom::tidy(lh_fit)
lh_metrics# Anova for Residual SD
lh_anova <- broom::tidy(anova(lh_fit))## Warning in tidy.anova(anova(lh_fit)): The following column names in ANOVA
## output were not recognized or transformed: Deviance, Resid..Df, Resid..Devlh_anova# Confidence Intervals
lh_conf <- confint(lh_fit)## Waiting for profiling to be done...lh_conf## 2.5 % 97.5 %
## (Intercept) -0.5209805 0.3235663
## Tarsus_mm 0.3599015 0.3859520# Coefficients & Credible Intervals
bayes_metrics <- broom::tidy(bayes_fit,intervals = TRUE)
bayes_metricsGeneral Observations from Exercise 5.3: 1) The confidence intervals, estimates, and std.error for all three fits are very close to one another. 2) If I didn’t initally clean the dataset (i.e. remove rows with nan), I would have proceeded with caution in interpreting the likelihood and bayesian fits. 3) I would also be hesistant of interpreting the bayes fit until I understood the prior assumptions.
5.3 Everyday I’m Profilin’
For your likelihood fit, are your profiles well behaved? For just the slope, use grid sampling to create a profile.
You’ll need to write functions for this, and use the results from your glm() fit to provide the reasonable bounds of what you should be profiling over (3SE should do).
Is it well behaved? Plot the profile and give the 80% and 95% CI. Verify your results with profileModel.
# Likelihood Function
likFun <- function(slope, intercept, resid_sd){
# Data Generating Process: Quail Fit for Varying Parameters
quail_fit <- intercept + slope * quail$Tarsus_mm
#Likelihood Function Based on a Log Normal Distribution
sum(dnorm(quail$Culmen_mm, quail_fit, resid_sd, log=TRUE))
}
# Determine Resid_SD bounds wtih 3SE from GLM Output
SE_bounds <- c(lh_metrics$estimate[2]-(3*lh_metrics$std.error[2]), #low bound
lh_metrics$estimate[2]+(3*lh_metrics$std.error[2])) # high bound
# Combinations for Likelihood Function
quail_lh_samp <- crossing(slope = seq(SE_bounds[1], SE_bounds[2],.001),
intercept = lh_metrics$estimate[1],
resid_sd = seq(0, 2, 0.01)) %>%
rowwise() %>%
mutate(logLik = likFun(slope, intercept, resid_sd)) %>%
ungroup()
# Likelihood
quail_slope_profile <- quail_lh_samp %>%
group_by(slope) %>%
filter(logLik == max(logLik)) %>%
ungroup()
# Quick Plot of Likelihood
qplot(slope,logLik, data = quail_slope_profile)
# Confidence Intervals - 95% and 80%
quail_slope_ci_95 <- quail_slope_profile %>%
filter(logLik > max(logLik) - (qchisq(0.95,1)/2)) %>%
filter(row_number()==1|row_number()==n())
quail_slope_ci_80 <- quail_slope_profile %>%
filter(logLik > max(logLik) - (qchisq(0.8,1)/2)) %>%
filter(row_number()==1|row_number()==n())
# Plot Confidence Intervals
lh_ci_plot <- ggplot(data = quail_slope_profile, aes(x = slope, y = logLik)) +
geom_line() +
geom_vline(data = quail_slope_ci_95, aes(xintercept = slope, color="blue")) + #80%
geom_vline(data = quail_slope_ci_80, aes(xintercept = slope, color="red")) + #95%
labs(color = "Confidence Intervals") +
scale_color_manual(labels = c("95%", "80%"), values = c("red", "blue")) +
xlim(SE_bounds) +
ylim(-1259,-1249) +
xlab("Slope") +
ylab("Log Likelihood")+
theme_bw()
lh_ci_plot## Warning: Removed 28 rows containing missing values (geom_path).
# Use Profile Model to Validate Function
lh_prof <- profileModel(lh_fit,
objective = "ordinaryDeviance")## Profiling for parameter (Intercept) ... Done
## Profiling for parameter Tarsus_mm ... Done# Plot Profile with Profile Model
plot(lh_prof, print.grid.points = TRUE)
# Look at CIs
lh_prof_ci_95 <- profileModel(lh_fit,
objective = "ordinaryDeviance",
quantile = qchisq(0.95, 1))## Preliminary iteration .. Done
##
## Profiling for parameter (Intercept) ... Done
## Profiling for parameter Tarsus_mm ... Doneplot(lh_prof_ci_95)
lh_prof_ci_80 <- profileModel(lh_fit,
objective = "ordinaryDeviance",
quantile = qchisq(0.8, 1))## Preliminary iteration .. Done
##
## Profiling for parameter (Intercept) ... Done
## Profiling for parameter Tarsus_mm ... Doneplot(lh_prof_ci_80)
confint(lh_fit)## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) -0.5209805 0.3235663
## Tarsus_mm 0.3599015 0.3859520Geneal Observations for Exercise 5.3 1) The CIs generated from the likelihood function were more constrained than the CIs generated from profileModel 2) ProfileModel CI results included the generated CIs and the differences were not alarming
5.4 The Power of the Prior
This data set is pretty big. After excluding NAs in the variables we’re interested in, it’s over 766 lines of data! Now, a lot of data can overhwelm a strong prior. But only to a point. Show first that there is enough data here that a prior for the slope with an estimate of 0.4 and a sd of 0.01 is overwhelmed by the data by demonstrating that it produces similar results to our already fit flat prior.
# Check priors used with bayes fit
prior_summary(bayes_fit)# Compare brm bayes to stan bayes computation
# Dinosaur Slow
bayes_fit_prior <- brm(Culmen_mm ~ Tarsus_mm,
family = gaussian(),
data = quail,
prior = prior(normal(0.4,0.01)),
file = "bayes_fit_prior.Rds")## Compiling the C++ model## In file included from file6122d57dc44.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Core:531:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122d57dc44.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:2:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/LU:47:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122d57dc44.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:3:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Cholesky:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Jacobi:29:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122d57dc44.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:3:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Cholesky:43:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122d57dc44.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/QR:17:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Householder:27:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122d57dc44.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:5:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/SVD:48:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122d57dc44.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Geometry:58:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122d57dc44.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:7:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Eigenvalues:58:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122d57dc44.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:96:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Sparse:26:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/SparseCore:66:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122d57dc44.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:96:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Sparse:27:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/OrderingMethods:71:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122d57dc44.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:96:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Sparse:29:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/SparseCholesky:43:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122d57dc44.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:96:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Sparse:32:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/SparseQR:35:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122d57dc44.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:96:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Sparse:33:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/IterativeLinearSolvers:46:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
## #pragma clang diagnostic pop
## ^
## In file included from file6122d57dc44.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:44:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp:14:13: warning: unused function 'set_zero_all_adjoints' [-Wunused-function]
## static void set_zero_all_adjoints() {
## ^
## In file included from file6122d57dc44.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:70:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/autocorrelation.hpp:18:8: warning: function 'fft_next_good_size' is not needed and will not be emitted [-Wunneeded-internal-declaration]
## size_t fft_next_good_size(size_t N) {
## ^
## 15 warnings generated.## Start sampling# Holy Cannoli Quick
bayes_fit_prior_stan <- stan_glm(Culmen_mm ~ Tarsus_mm,
family = gaussian(),
data = quail,
prior = normal(0.4,0.01))
# Compare Bayes Fit Priors Against Bayes Fit
posterior_interval(bayes_fit)## 2.5% 97.5%
## b_Intercept -0.5050770 0.3020018
## b_Tarsus_mm 0.3602278 0.3855762
## sigma 1.1794704 1.3022714
## lp__ -1260.0431710 -1255.4884869posterior_interval(bayes_fit_prior, digits=5)## 2.5% 97.5%
## b_Intercept -0.7257101 -8.511966e-03
## b_Tarsus_mm 0.3703985 3.921486e-01
## sigma 1.1798552 1.306030e+00
## lp__ -1258.8974800 -1.254359e+03posterior_interval(bayes_fit_prior_stan, digits=5) ## 5% 95%
## (Intercept) -0.9212210 -0.5150667
## Tarsus_mm 0.3865909 0.3982382
## sigma 1.1947696 1.2956723Second, see if a very small sample size (n = 10) would at least include 0.4 in it’s 95% Credible Interval. Last, demonstrate at what sample size that 95% CL first begins to include 0.4 when we have a strong prior. How much data do we really need to overcome our prior belief? Note, it takes a long time to fit these models, so, try a strategy of spacing out the 3-4 sample sizes, and then zoom in on an interesting region.
sampFun <- function(samp_size = 10, predictor = quail$Tarsus_mm, response = quail$Culmen_mm, prior_slope = 0.4, prior_slope_sd = 0.01){
# Make Data Frame of Predictor and Response Variables
df <- data.frame(predictor, response)
# Extract Predictor and Response Variable
sampledData <- df[sample(nrow(df), size = samp_size, replace = FALSE),]
#Likelihood Function Based on a Log Normal Distribution
bayes_fit_prior_stan <- stan_glm(response ~ predictor,
family = gaussian(),
data = sampledData,
prior = normal(prior_slope,prior_slope_sd))
# Gather CI Information
bayes_metrics <- posterior_interval(bayes_fit_prior_stan)
return(bayes_metrics[2,])
}
# Test out Function for a Small Sample Size (N = 10)
sampFun(10, quail$Tarsus_mm, quail$Culmen_mm)## 5% 95%
## 0.3925685 0.4071586# Test Out Function on Different Sample Sizes
# Sample Size = 100
sampFun(100)## 5% 95%
## 0.3921760 0.4055417# Sample Size = 200
sampFun(200)## 5% 95%
## 0.3914971 0.4045587# Sample Size = 300
sampFun(300)## 5% 95%
## 0.3904886 0.4036722# Sample Size = 500
sampFun(500) # BINGO## 5% 95%
## 0.3876449 0.3999444General Observations for Example 5.4 1) Bayes fits with large sample sizes overwhelm the prior 2) Bayes fits with small sample sizes include the prior (0.4) 3) A sample size of 500 through sampling without replacement dipped belows 0.4 4) A more efficient way to write a function would have been to determine from the CI if it contained CI the prior and if it contained the prior, it would skip to a larger sample size until the prior was no longer included in the CI. I’m about half way there, but I figured I should turn this in since it’s getting close to midnight.
6. Extra Credit
I obviously didn’t do this one, but I have a question: Could you use Bayesian stats to determine if a close election was going to a recount or runoff? For example, Stacey Abrams vs. Brian Kemp. Stacey Abrams has not conceded because she is waiting for every vote to be counted, which includes absentee votes and all counties voting. With each day more votes are being reported; thus, I question if you can use Bayesian stats to 1) Estimate the vote count 2) determine how close vote are between candidates.